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Abstract 

The well known linear relation between density and peculiar velocity distributions is a powerful 
tool for studying the large-scale structure in the Universe. Potentially it can test the gravitational 
instability theory and measure Q. At present it is used in both ways: the velocity is reconstructed 
provided the density is given and vice versa. Reconstructing the density from the velocity field 
usually makes use of the Zel’dovich approximation. However, the standard linear approximation in 
Eulerian space is used when the velocity is reconstructed from the density distribution. I show that 
the linearized Zel’dovich approximation, in other words the linear approximation in the Lagrangian 
space, is more accurate for reconstructing velocity. In principle, a simple iteration technique can 
recover both the density and velocity distributions in Lagrangian space, but its practical application 
may need an additional study. 


1 Introduction 


In this talk I would like to discuss the quasi-linear regime of the gravitational instability and in 
particular the density-velocity relation. This is the relation between peculiar velocity v p (x) and 
density contrast <$(x). The peculiar velocity is defined as 


Vp(x,f) 



( 1 ) 


where x = r/a is a comoving coordinate and a = a(t ) is the scale factor. The density contrast is the 
relative density perturbation 


(2) 


The equation relating peculiar velocity v p and the density contrast 6 is derived, for example, in [8]) 



( 3 ) 


where H = a/ a is the Hubble constant, / = and D is the growth factor of density fluctuations 

6 a D(t)\ in the ft = 1 matter dominated cosmology D oc a oc < 2 / 3 . Eq.(3) is widely used for 
reconstruction of velocity field from density field (see e.g., [10]). Quantitatively eq.(3) is valid only 
in the linear regime which is based on two assumptions: (1) 6 << 1 and (2)x ~ q, where q is the 
unperturbed (Lagrangian) comoving coordinate of the particle. The generalization of relation (3) 
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including the nonlinear corrections to b was studied in [7]. The second assumption is a subject of the 
paper. To stress this assumption we obtain an analog of eq.(3) from the Zel’dovich approximation 
[12], [11]. The Zel’dovich approximation is usually formulated as a relation between the comoving 
Eulerian x and Lagrangian q coordinates of particles 

x = q + D(t) ■ s(q) , (4) 

where s(q) = - V V’o(q) * s the potential vector field characterizing the growing mode of the initial 
perturbations. The density can be found in terms of the eigenvalues Ai(q), A 2 (q), A 3 (q) of the initial 
deformation tensor field d{ t(q) = —dsi/dqk 

piq,t) = (1 - £A0(T - £>A 2 )(1 - da 3 ) (5) 

and the peculiar velocity field in terms of s(q) 

v p (q,t) = aDs(q). (6) 


In the limit | D A, |<< 1 eq.(5) reduces to 




( 7 ) 


where dsi/dqi = Ai + A 2 + A 3 . Solving eq.(7) together with eq.(6) we readily obtain the desired 
Lagrangian analog of eq.(3) 


v p (q,o) = 


H fa 
47T 



q' - q 
|q'-q | 3 



(8) 


The only difference between (3) and (8) is the difference between Eulerian coordinates x and La- 
grangian coordinates q. One might think that at small 6 the difference between Eulerian and La- 
grangian coordinates always can be neglected. This assumption being right for some spectra is gener- 
ally wrong and it is the major subject of this talk. 


2 Three regimes of gravitational instability 


A very useful theoretical model for studying the formation and properties of the large scale structure is 
a dust like continuous medium described in great detail in [8]. The equations describing the evolution 
of density perturbations in this model can be conveniently written in a slightly different form 


dy diy-Vi) _ 0 

dD + dxi 


dvi dvi 
dD + Vk 8x ~ k 


3 tl 


2Dp y dxi 


. dtp 

(^r- + u *) 


d 2 q> _ s 

dx 2 x ~ D 


V = fj( 1 + 6), 


( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 


where as usually x, = r;/a is the comoving coordinate, but the peculiar velocity is scaled by the 
first derivative of the growth factor v t = u pt /(a • D) — D~ x • dxi/dt = dxi/dD, y = p • a 3 , and the 
perturbation of the gravitational potential <p g = 3/2 Qa 2 D ■ <p. 
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The standard gravitational instability paradigm usually assumes that the evolution of the density 
perturbations having a comoving scale L has two stages: the linear stage when cr p (L) < 1 or equiv- 
alently L > 1 / k n i and afterward the nonlinear stage when a p {L) > 1 and L < 1 /k n i, where k n [ is 
defined by the condition 

kni 

a p(k~l) = 4tt • D 2 J P(k)k 2 dk = 1. (13) 

o 

I show that for some initial spectra an additional stage which can be called the quasi-linear regime 
can be identified. For such a spectra the difference between the Eulerian (r t ) and Lagrangian (<?,) 
coordinates is significant and can influence the density-velocity relation. 


2.1 Linear regime 


The evolution of the scales which are in the linear regime can be described by a set of simple linear 
equations. Linearizing eq.(9) and (10) we obtain 


86 

8D 



^o(x), 


(14) 



The well known growing solution of this system in our variables takes a form 


(15) 


6(x,D) = D-6o{x), (16) 

r,(x,Z?) = u 0 .(x), (17) 

¥>(x) = ¥>o(x), (18) 

where 6 0 (x), v 0 ,(x) and y?o(x) are the initial density, velocity and gravitational potential (scaled as 
was indicated above) perturbations. The growing mode is specified by only one spatial function, so 
the initial density, velocity and potential are related as voi(x) = — and 6o(x) = — From this 
solution one can see that it preserves the initial spatial structure of the perturbation. In other words 
this approximation ignores the fact that the growth of the amplitude of the perturbation requires the 
displacement of mass. The linear approximation can be used quantitatively only when 6 <C 1, but in 
practice it is often pushed to the limit 6 « 1. However, a simple order of magnitude analysis suggests 
that in some cases there is an additional condition restricting the linear approximation in the above 
form. 

The continuity equation (eq.(9)) can also be written in terms of the density contrast 


86 8v, 86 8v{ 

8D dxi V ' 8x{ + dxi 


0 . 


(19) 


For the order of magnitude estimates we replace the derivatives by the finite ratios 


6/D + v/l v + v(6/l s ) + 6(v/l v ) ~ 0. 


( 20 ) 


Here v is the characteristic fluid velocity and l$,l v are the coherence lengths for spatial variations of 6 
and v respectively. Comparing linear eq.(14) with the exact continuity equation (19) one can see that 
two last terms were discarded from eq.(19). The order of magnitude estimate of eq.(14) simply yields 
6 ~ D • ( v/l v ) which is of course in agreement with solution (16). 

Making use of eq.(20) we readily find that the last term in eq.(19) can be neglected if the density 
contrast is small 6 < 1 or equivalently D < l v /v. 
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However, discarding the third term of eq.(19) requires a condition 6 < ls/l v or equivalently D < l$/v 
which may be much stronger: if Is <C l v then the linear theory is correct only until 6 < ls/l v C 1. 

Discarding the second term on the left hand side of the Euler equation (eq.(10)) is justified if 
D < l v /v which coincides with the requirement for the density contrast to be small. The right 
hand side of the Euler equation is zero to the linear order. Of course, using the linear theory for a 
quantitative analysis all signs “<” in the above inequalities must be replaced by 

In a simple case when Is ~ l v the linear approximation can be used for a qualitative analysis up 
to the beginning of the nonlinear stage 6 « 1. The relation between Is and l v is determined by the 
initial spectrum and will be discussed later. 

The Fourier components <5 k of the perturbations in the linear regime grow linearly and the phases 
remain constant 

£ k (D) « D • <^(0) (21) 

tM-D) ~ V’k(O)- (22) 


2.2 Nonlinear regime 


The perturbations with scales smaller than k n * are in th e nonlinear regime. The analysis of the 
nonlinear perturbations requires solving the full nonlinear equations (9-12) which at present is not 
possible in a general case. In a few cases when a particular type of symmetry is imposed (e.g. 
spherical) the exact solutions are known. For the random initial conditions the hierarchical clustering 
model gives a good qualitative and in some cases quantitative description of the process (for details 
see, e.g. [8] and references therein). In the nonlinear regime both equalities eq.(21) and (22) break 

6 k (D)?D-6 k ( 0) 

^k(£)#tM0). 

The growth of the amplitude slows down and the growth of the density perturbations is determined 
mostly by the phase adjusting. 


2.3 Quasi-linear regime 

Two nonlinear terms dropped from the continuity equation (eq.(19)) and Vk§^ dropped from 

the Euler equation (eq.(10)), when linearizing the original equations, can be relatively easy retained 


86 86 

8D + V '8xi 


6 o(x) 


8v , dvi _ 

dD +Vk 8^-°- 

(^(x) = <y?o(x). 


This set of partly linearized equations differs from one usually used in the linear analysis. In the 
Lagrangian form these nonlinear terms can be absorbed into the time derivative retaining linearity 


7b " ‘“O' 


The solution to this system is obvious: 



6(x,D) = D-6 0 ( q), 


(23) 
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Vi(x,D) = t? w (q). (24) 

Solution (23, 24) looks similar to (16, 17), but must be interpreted differently. It does not describe the 
’real’ distribution of density and velocity in space, instead it describes them in the Lagrangian space. 
However, we may find the ’real’ (Eulerian) distributions assigning the density and velocity values of 
eq.(23, 24) to the coordinates x, = qi + D ■ u 0 i(q), provided that x, = g, at D = 0. The last equation is, 
of course, the Zel’dovich approximation eq.(4) with s, = vo». The difference between the Lagrangian 
coordinate q t and the Eulerian coordinate X{ can be unimportant, for instance, if we are interested 
in the internal structure of clumps, or their integral properties like masses, or angular momenta. 
However, it is crucial if the goal is the spatial distribution of the clumps or the velocity on scales 
of tens of Mpc.. Passing by we note that the Zel’dovich approximation is actually an extrapolation 
of solution (4) until the shell crossing. Formally it is not accurate even to the second order, but 
practically is quite good [2], [11]. 

The difference between Eulerian and Lagrangian coordinates comprises the difference between the 
linear and quasi-linear regimes. 

In terms of the Fourier components the qasi-linear regime is a mixture of the linear and nonlinear 
cases. The Fourier components grow approximately according to the linear theory (eq.(21)) and the 
phases are very different from the initial ones 

6 k (D)*D-6 k (0), 

MD) # 1&k(0). 

An immediate question arises: which scales are in the quasi-linear regime? The answer comes 
from our order of magnitude analysis: the scales in the quasi-linear regime are those which on the 
one hand are greater than the nonlinear scale eq.(13) and on the other hand are smaller than the 
typical displacement of a particle from its unperturbed position (Lagrangian coordinate) to the current 
position (Eulerian coordinate). Ryden and Gramann [9] found in two-dimensional N-body simulations 
that the scale where the phases become substantially different from the initial ones in some cases 
does not scale with k n i. As the author checked it perfectly scales with the theoretically calculated 
characteristic displacement of mass. 


3 Characteristic displacement of mass 

For estimation ot the mean distance passed by the mass one can use the Zel’dovich approximation 
eq.(4) 

d rm , =< (x - q) 2 > 1/2 = D(t) ■ s rms 

vj 1ipt*p 

Srms =< s 2 (q) > 1/2 (25) 

However eq.(4) holds only before crossing of the orbits. To make use of relation (4) at later times one 
can smooth initial perturbation field at k n i which is defined by eq.(13). Then 

knl 

s 2 rms (k nt ) = 4* J P(k)dk. (26) 

0 


Combining eq.(25), (13) and (26) we may express the r.m.s. distance in terms of the initial power 
spectrum 


,(knl) 


P(k)k 2 dk 


(27) 
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Obviously integrals in eq.(27) could be easily defined in terms of an arbitrary (e.g., Gaussian) smooth- 
ing window function. Pushing the upper limits in eq.(27) to the infinity one obtains the familiar 
characteristic scale of the initial vector field s(q). On the other hand, the r.m.s. displacement is 
always equal to the characteristic scale of the initial velocity field smoothed with the scale on nonlin- 
earity. Before going to implications of eq.(27) it is worth discussing its accuracy. One obvious test of 
eq.(27) is its comparison with N-body simulations. This test has been carried out for a series of power 
law models. 


3.1 Power law spectra 


Pure power law initial spectra play a very important theoretical role since they have various scaling 
properties. Here I discuss slightly modified power law models imposing explicitly a cutoff k i at small 
k corresponding to the large scales. This modification solves the problem of divergence at large scales 
for certain spectra. The cutoff is arbitrary and can be treated as a free parameter. The limiting case 
of a pure power law can be analyzed assuming an obvious limit k\ — ♦ 0. All iV-body simulations 
assume a sharp cutoff of the initial spectrum at the fundamental mode kf = 2n/Lb 0X . Therefore this 
theoretical model better suits the purpose of explaining the results of ,/V-body simulations. 

The initial spectrum has a form 



k > k\ 
k < k i 


(28) 


where A n is a constant. Assuming the linear growth of the amplitude of density perturbation at 
k < k n i eq.(21) one can easily find the nonlinear scale 




Assuming k n i >> k\ one obtains a familiar result [8] 

*ni «zr 2/(n+3) ;rc> -3, 


which tells us that the scale of nonlinearity monotonically grows with time if n > -3. This is usually 
interpreted as hierarchical clustering of smaller nonlinear clumps into larger ones. The typical mass 
of clumps formed can be estimated as 

M nt ~ p^-k-f (X T) 6 /("+ 3 ); n > -3. (29) 

It is also assumed that a substantial fraction of mass (~ 0.5) is in such clumps. Obviously, for making 
such a clumps it is sufficient to move mass over a distance about k~f . On the other hand applying 
eq.(27) with the lower limits at k = Aq to spectra eq.(28) one easily obtains in the limits k n i >> k x 


t 




fn ± 3 1 

V w+1 k nl 
J2tn- 


1 

knl ’ 




71 > — 1 

n = — 1 


-1 > n > — 3 


(30) 


At n > -1 d Tma is about k n i, however at n < -1 the r.m.s. distance can be much greater than k ~ * 
depending on k\. This means that in the case -3 < n < -1 the clumps with masses ~ M n i (eq.(29)) 
move coherently, and as we commented before the coherence scale of this motion always about d rma . 
If the initial spectrum is even steeper n < — 3 then no clumps form, until the first pancakes are formed. 
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3.2 Comparison with 7V-body simulations 


Prediction of eq.(27) was tested in a series of 3 dimensional TV -body simulations. The simulations 
[3] were done with 128 3 particles on equal mesh for five power law initial spectra eq.(28) with n = 
1,0, — 1, —2, —3 and k\ = kj = 2ir/Lb ox each model was simulated with three different sets of random 
numbers. Each simulation was stopped when k n \ = 64,32,16,8 and 4 and d rms (N — body ) was 
calculated from the particle distribution 

drms (n - body) =< (r, - q ,) 2 > 1/2 . 

The results in mesh units are shown in Figure 1. Along with the experimental points the analytic calcu- 
lations using eq.(27) with the lower limit k = kj are also shown with dashed lines (n = 1,0, — 1, -2, -3 
from the bottom to top). The solid line shows relation d rms = Three symbols of the same type 
at given k n \ corresponds to different realization of random numbers. For n = 1 (open triangles) and 
n = 0 (open squares) the results for different realizations almost coincide with each other so fewer than 
three symbols can be seen. The n = —3, —2, —1 series are almost exactly coincide with the theoretical 
predictions for all k n \ from 64£/ to 4 kj. The agreement between the theory and the experiments gets 
worse for the n — 0 and especially for the n = 1 series. Actually it is not unexpected. According to 
eq.(30) for the n = 1 and n = 0 models d rms is only 1.4 and 1.7 times greater than therefore the 
accuracy of the essentially linear theory can not be very good. However, even for these spectra in the 
worst case the discrepancy does not exceed about 40%. Applying the theory to cosmology one has 
to keep in mind that all realistic models of the large scale structure based on gravitational instability 
suggest that the present nonlinear scale is in the region where the initial spectrum had effective slope 
about n ~ — 1 or even steeper. For such spectra the accuracy of eq.(27) is far better. In any case it 
gives a good lower limit for d rms . 

3.3 CDM spectrum 

The initial spectrum of perturbations might not be a pure power law. A good realistic test model is 
the CDM spectrum. In this case a numerical evaluation of the integrals in eq.(27) is needed. The 
typical normalization k~l » 4.5 h~ l Mpc roughly corresponding to the COBE normalization of the 
primordial spectrum gives d rms (C DM) % 9 h~ l Mpc & 900 km/s. 


4 Discussion 

The growth of the density perturbations necessarily assumes the displacement of mass. A good 
estimate of the typical distance traveled by a particle (e.g. a galaxy) from its unperturbed position can 
be given analytically applying the Zel’dovich approximation to the initial power spectrum truncated 
at the nonlinear scale eq.(27). The comparison with the three-dimensional N-body simulations carried 
out for five power law models with n = —3, —2, — 1,0, 1 showed that eq.(27) is almost exact if n < — 1 
and for n = 0 and 1 the actual displacement is about 20-40% greater than the analytical prediction. 
Actually, the density distribution predicted by the truncated Zel’dovich approximation is far better 
than all other tested analytical or semianalytical approximation (see [1], [4], [5], [6]). 

For the CDM spectrum (taken here as an example) with the COBE normalization the theory 
predicts the characteristic displacement of about 900A;m/s. The statistic of the displacement is roughly 
Gaussian, which means that more than 30% of the mass as well as of galaxies has been displaced even 
more. It implys that the standard linear relation between 6 and v p eq.(3) can result in a large error 
unless the smoothing scale is greater than at least lOOOfcm/s. 

Smoothing effectively reduces the difference between the Eulerian and Lagrangian coordinates and 
with a sufficiently large smoothing scales it can be completely erased. However, the displacement 
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pattern due to gravitational instability is not identical to that of smoothing, therefore a reasonable 
quantitative agreement between the Eulerian and Lagrangian distributions may require even a larger 
smoothing scale. For example, the analysis of the nonlinear Eulerian and correspondingly scaled initial 
Lagrangian density distributions smoothed with the same Gaussian window showed that only when 
the size of the window reaches R w « 1.65d rmi the r.m.s. difference between the two density fields 
becomes about 25% almost independently of the initial spectrum [3]. Applying this result to the 
CDM model we find that for obtaining a density field for the density- velocity relation eq.(3) with 
25% accuracy the original density field must be smoothed with a scale at least 1500fcm./ s. 

A similar problem arises when a reconstruction technique is tested in an N-body simulations. If the 
r.m.s. displacement of particles in the simulation is smaller than in the tested theory or the universe 
(e.g. because the box is too small) then the technique will look better in the numerical test compared 
to its performance in the theory or on real data. 

The above constraint seems not to be a problem for the POTENT which employs the Zel’dovich 
(or modified ZePdovich) approximation, but for the reconstruction of the peculiar velocity field from 
the density distribution may suffer from this effect. 

An obvious way of eliminating this effect consists in recovering the Lagrangian coordinates. One 
may start from the original density field and find the first approximation to the peculiar velocity field 
using eq.(3). Next, the density must be assigned to the Lagrangian coordinates which can be found 
from the reverse ZePdovich approximation g, = X{ - D • v t . After that, the whole procedure must be 
repeated again, but this time eq.(8) must be used. Iterating, one can get the Lagrangian distributions 
with desired accuracy. Unfortunately this method requires the knowledge of the factor D which in 
turn depends on f l and the normalization of the primordial spectrum. 
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Fig.l d rms as a function of the nonlinear scale. 

Points are results of the N— body simulation. 
Dashed lines are the theoretical curves. 

A solid line shows d^, = 


l/k nl . 




